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ABSTRACT 



The Pseudo Wigner-Ville Distribution is a signal transformation of an input time signal 
into a joint time-frequency domain, that provides an excellent characterization of an input 
signal as well as it’s respective energy content. A smoothing window and energy sensitivity 
analysis of the discretized Pseudo Wigner-Ville Distribution is presented along with a 
machinery monitoring application using the resulting Wigner-Ville energy of the sampled 
signal. The ability to apply the Pseudo Wigner-Ville Distribution to both steady-state and 
transient machinery may enable the monitoring and diagnostics of virtually any piece of 



equipment. 
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I. INTRODUCTION 



The recent world events have caused drastic changes in many areas and made the future 
that much more unpredictable. With this high degree of uncertainty and drastic financial 
cutbacks as evidenced in the Defense budget, the emphasis of technological advance and 
dominance is all the more heightened. With the requirement to "do more with less", it is 
imperative to find ways to save money. It is with this mind-set that better ways to prolong 
the life of operating machinery and defray unneeded maintenance costs be investigated. For 
many years, most transient machinery has gone unmonitored and removed at periodic time 
intervals for refurbishment regardless of the current machinery condition. Developing a 
method to monitor and diagnose these machines of a transient nature will provide a way to 
determine when overhaul or replacement is needed. 

Vibrational data has long been used as the tool to measure a machine’s health. Early 
methods of monitoring machinery had been for the experienced personnel to listen or touch 
the machine to detect if a fault exists, however, this method is dependent on human senses 
and the same individual conducting the monitoring. It was for this reason that vibration levels 
as monitored by a dynamic signal-analyzer are now used for machinery health monitoring. 
The vibration levels of a good condition establish a baseline to allow a means of comparison 
when subsequent measurements are taken. These dynamic signal analyzers are time 
independent and do not show the time dependency required by transient machinery. It is 
because of this limitation, that a three-dimensional time, frequency and energy representation 
such as the Pseudo Wigner-Ville Distribution is being investigated as a method to monitor 
transient machinery. 
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A. PREVIOUS WORK AND ACCOMPLISHMENTS 



This thesis is a follow-on work which has been preceded by LT G. Rossano and LCDR 
J. Calahorrano, both in conjunction with Prof. Y.S. Shin. Rossano established the initial form 
of the computer program that calculates the discretized Pseudo Wigner-Ville Distribution. 
He evaluated the effect that different signals and procedures had on the Distribution along 
with characterizing an actual pump signal. Calahorrano made some further improvements to 
the computer program and then applied this improved version of the computer program to 
characterize some classified signals. Both of these works have been a great benefit in dealing 
with characterizing signals, but have not investigated the smoothing method and energy or 
amplitude of the discretized Pseudo Wigner-Ville Distribution. 

B. OBJECTIVES OF CURRENT RESEARCH 

The Pseudo Wigner-Ville Distribution provides an outstanding time-frequency domain 
spectra representation of an input signal as shown through the works among others of Rossano 
and Calahorrano here at the Naval Postgraduate School. The purpose of this research is: 

• To review and modify the computer program developed here at the Naval Postgraduate 
School to support the smoothing window and energy analysis work. 

• To conduct a smoothing window sensitivity analysis. 

• To investigate the energy content of the discretized Pseudo Wigner-Ville Distribution. 

• To apply the developed smoothing window and energy principles to an actual 
machinery monitoring application. 

• To investigate the feasibility of using the Pseudo Wigner-Ville Distribution energy as 
a pre-processed input to a Neural Network. 
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II. THE PSEUDO WIGNER-VILLE DISTRIBUTION 



A. THEORETICAL BACKGROUND 
1. Evolution 

The current day Pseudo Wigner-Ville Distribution is a three dimensional (time, 
frequency, and energy) representation of a signal that is particularly well suited for analysis 
of non-stationary signals. Eugene Wigner [Ref. 1] first introduced the Wigner distribution in 
1932 to study the problem of statistical equilibrium in the area of quantum mechanics. His 
work was furthered in 1948 by J. Ville [Ref. 2], who used the Wigner distribution in the area 
of signal analysis. A major contribution of Ville’s, is the use of analytic signals as an input 
to the distribution vice the customary real signal. The advantages of using an analytic signal 
in the distribution are two-fold [Ref. 3]. First, the distribution of a real signal results in a 
symmetrical spectrum with only half of the result containing useful information, while the 
use of an analytic signal avoids this negative frequency redundancy. Secondly, by accounting 
for only the positive frequencies it satisfies both the practical and the mathematical 
completeness of the problem. The third part of the Pseudo Wigner-Ville Distribution (ie. 
Pseudo) , arises from the need to apply a smoothing window to the resulting distribution. 
This is a result of the presence of cross terms that arise from multiple component signals and 
will be covered in depth later. 

The works of T.A.C.M. Claasen and W.F.G. Mecklenbrauker [Refs. 4,5, and 6] in 
1980 have paved the way for many recent applications of the Pseudo Wigner-Ville 
Distribution. Their three part paper is an all encompassing work that has served to highlight 
the capabilities of the Pseudo Wigner-Ville Distribution and allow for greater knowledge of 
the subject with an emphasis on signal processing and analysis. 
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A considerable amount of recent work using the Pseudo Wigner-Ville Distribution 
have come in the areas of optics [Refs. 7,8, and 9] and speech [Ref. 10 and 1 1]. Additionally, 
T.J. Wahl and J.S. Bolton of Purdue University [Ref. 12] have used the Pseudo Wigner-Ville 
Distribution to analyze structural impulse responses. Two recent papers have investigated the 
use of the Wigner-Ville Distribution in the area of machinery monitoring. Flandrin et. al. 
[Ref. 13] proposed the Wigner-Ville as a means to confirm a machinery diagnosis in the 
specific application of a Pressurized Water Reactor power plant and lastly, Forrester [Ref. 14] 
has investigated the Wigner-Ville Distribution as a method for fault detection in gears. As 
evidenced by the cited examples, the capability and adaptability of the Wigner-Ville 
Distribution is enormous with applications of it’s use rapidly expanding. 

2. Function Definition 

The Wigner-Ville Distribution is a transformation of a signal into the time- 
frequency domain. By definition, the cross Wigner-Ville Distribution is defined as: 

o© 

WDF Z ' S - fe-j"' r(t+-|0 s*(t — dx 

J £ £ 

— O© 

where: r(t) = a complex time history 

s(t) = a complex time history 
t = time 
w = frequency 
* = complex conjugate 

Similarly, the auto Wigner-Ville Distribution is defined as: 

oo 

WDF Z fe- Jax r(t+4> r*(fc-— ) dx 

r ' r J 2 2 
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Since the concentration of this research is with the representation of a single input signal, the 
auto Wigner-Ville Distribution will be used. This expression is essentially a Fourier transform 
of the auto-correlation of a signal, which may be thought of as a three dimensional spectral 
density function. 

3. Distribution Properties 

There are several properties of the Distribution that are worth noting. The 
properties listed below have been shown in many works concerning the Wigner-Ville 
Distribution. These properties will be shown in the form of the auto Wigner-Ville 
Distribution [Ref. 15]. 

Time and frequency shift properties follow: 

• A time shift in the signal corresponds to a time shift of the Distribution: 

WDF s(t^)sU-^) (t ' “ WDF SiS U~ T,W) 

• A frequency shift in the signal corresponds to a frequency shift in the Distribution: 

WDF el az s e ]a Cs ( t, 0)) - WDF s s (t, U-Q) 

• Both a frequency and time shift in the signal corresponds to the same in the 
Distribution: 



(t, «} - WDF S ' S (t-T, W-Q) 



Since our concern is in the area of machinery monitoring and diagnostics, the aforementioned 

nece^f'i y 

properties are a nessesity in the development of a suitable monitoring program. 
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The next three properties provide information concerning the energy contained 
within the Wigner-Ville Distribution. 

• The integration of the Distribution over the frequency domain provides the 
instantaneous signal power: 



- f WDF S _(t,o>) dw - I sit) I 2 
2n J s,s 



• The integration of the Distribution over the time domain provides the energy density 
spectrum: 



f WDF s s (t,u) dt - I 5(o) Y 



♦ The integration over both the time and frequency domain provides the total energy of 
the input signal: 



— f f WDF . ( t, o ) dtdo - I I s I I 2 

271 J J 

— oo- oo 

Again, as with the time and frequency properties, the energy contained in the Wigner-Ville 
Distribution must be definable in order to use this in a machinery monitoring program. 
Additional properties may be found in Reference 4. 

4 . Cross Terms Within The Distribution 

A particularly troublesome characteristic of the Wigner-Ville Distribution is the 
presence of cross terms. These cross terms result in the Distribution taking on negative values 
as shown later in Figure 2. Since the properties above have shown that the Wigner-Ville 
Distribution represents an energy value, these negative values raise an area of concern. Jones 
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and Parks [Ref. 16] show that these cross terms are essentially the interference of multiple 
component signals that occur when there is a non-zero cross-correlation within the signal. 
The equation listed below shows the cross-term that results from a signal that consists of two 
components. 

WDF x+y ( t , 0) ) - WDF x (t, 0)) + 2RE[WDF xy (t,<*>)] + WDF y (t,u>) 



As shown, the cross-term is represented by the center term while the other terns represent 
their respective signal components. There are several methods for eliminating these terms and 
providing a more presentable Distribution and there has been a considerable amount of 
research in the are of an optimum smoothing method for the Wigner-Ville Distribution as 
indicated by References 3,4,17, and 18. The method chosen for smoothing the Wigner-Ville 
Distribution in this work is a post processing, sliding Gaussian window function. What is 
meant by post processing, is that the distribution is smoothed after calculation of the Wigner- 
Ville Distribution. 

5. Gaussian Smoothing Function 

The Gaussian window function was first proposed mathematically as a means to 
smooth the Wigner-Ville Distribution by N.D. Cartwright [Ref. 19]. The continuous form of 
the Gaussian window function follows: 



<?( t, o>) 




2 no t o M 



Wahl and Bolton in Reference 12 have furtherd this work and developed the sampled form 
of the window function that is used in our discretized form of the Pseudo Wigner-Ville 
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Distribution for smoothing. The effect of various size smoothing windows will be presented 
in greater depth in Chapter III. 

B. WIGNER-VILLE DISTRIBUTION COMPUTER PROGRAM 

The computer program listed in the Appendix is a method that was developed here at 
the Naval Postgraduate School for calculating the discretized Pseudo Wigner-Ville 
Distribution. As mentioned earlier, the program that computes the discretized Pseudo 
Wigner-Ville Distribution is largely the work of G. Rossano [Ref. 15], who initially 
investigated the use of the Wigner-Ville Distribution as a tool for machinery monitoring. For 
the purposes of this work, the program has been thoroughly reviewed and rewritten to allow 
for ease of readibility and modification. The current version enables the user to look at not 
only the resulting smoothed version of the Pseudo Distribution, but also the unsmoothed 
Wigner-Ville Distribution. The second major modification was to program the smoothing 
window subroutine to allow for the use of different size windows in smoothing the 
Distribution. These changes have proven invaluable for evaluating the effect of various sized 
smoothing windows and in conducting the energy sensitivity analysis. Presented below are 
the two major equations that comprise the calculation of the discretized Pseudo Wigner-Ville 
Distribution as performed in the program listed in the Appendix. 

1. Discrete Time Wigner-Ville Distribution 

As mentioned earlier, the Wigner-Ville Distribution is essentially a Fourier 
transform of an auto-correlation function. There have been a couple of equations developed 
to represent the discrete Wigner-Ville Distribution. The version that has been used in the 
development of the program listed in the Appendix, is one that was developed in Reference 
4 by Classen and Mecklenbrauker and and later used in the work of Wahl and Bolton [Ref. 
12] and Rossano [Ref. 15]. The only difference of note between the various discrete versions 
is essentially a difference in constants. All of the versions are basically a Fourier transform 



8 



of an auto-correlation function. Shown below is the discretized equation that has been used 
for programming purposes in this work: 



N 

WDF g s (iA< 0 , jAt) - £ ReiFFT[2A t CORR(i)]) 

j - 1 



where: Re = Real part 

Corr(i) = complex auto-correlation of the signal s(t) 
FFT = Fast Fourier Transform 
At = sample time 



2. Discretized Gaussian Window Function 

Due to the need to deemphasis the resulting cross-terms within the Wigner-Ville 
Distribution, a suitable method for smoothing has been obtained. The sampled equation of 
the selected smoothing window is shown below [Ref. 12]: 



where: 



G(t,w) 

-2M<m<2M & 

a t = MAt & ct,, 

lr W 



1 

2n ° t a <0 

-2N<n<2N 
= NAw 



-l 



e 



(mAO 2 

2 o 2 c 



(nAcj) 2 , 
2o 2 
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III. SMOOTHING WINDOW ANALYSIS 



As shown above, the smoothing function or window that is used in this work is a 
Gaussian window that is applied to the unsmoothed Wigner-Ville Distribution. Figure 1 is the 
raw time input signal of a 100 & 400 Hz sine wave that will be processed through the Wigner- 
Ville Distribution. Figure 2 is the unsmoothed Distribution that results from this combined 
100 and 400 Hz sine wave. As can easily be seen in Figure 2, the resulting cross terms are 
very dominant and need to be deemphasized. This Gaussian smoothing function is applied 
after the Wigner-Ville Distribution has been calculated (ie. post processing), which allows for 
visual inspection of both the unsmoothed Wigner-Ville Distribution and the effect that 
different sized windows may produce on the smoothed Distribution (ie. Pseudo). For 
example. Figure 3 is the smoothed result of Figure 2 after applying a 13 by 13 sized window. 
The next few sections will discuss some of the characteristics of this post-processing 
smoothing window and what may be gained or lost through varying the size of the applied 
window. 



10 



100 0 400 HZ SINE WAVES 
SIGNAL LENGTH = 0.256 SEC / AMPS = 1.0 




I : I I I r 

0'2 01 00 OV 02- 

(suoa) ganiiidiAjv 



Figure 1. 100 & 400 Hz Sine Wave Input 



11 



0.000 0.064 0.128 0.192 0.256 

TIME (sec) 



100 & 400 Hz SINE WAVE 
UNSMOOTHED DISTRIBUTION 
AMPLITUDES = 1.0 



- n 




Figure 2. Unsmoothed Wigner-Ville Distribution 
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igure 3. Pseudo Wigner-Ville Distribution (13 X 13 Window) 
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A. SIZE CRITERION 



As shown by N.D. Cartwright [Ref. 19] for the continuous form of the Wigner-Ville 
Distribution, the size of the smoothing window may be selected such that a positive 
Distribution will result. Following, are these criterion that have been mathematically proven 
for the continuous case to assure a positive distribution: 




where: 

a r - MA t a,, - i\7Aa> 

C to 

By chosing M and N such that this criterion is met, the resulting smoothed discretized 
Wigner-Ville Distribution does not necessarily result in all of the values being positive. This 
may be seen in Figure 4, where the window that was used was again a 13 by 13 (ie. N = 13 
& M = 13) sized window. However, it does deemphasize the large cross-terms to the point 
that they are nearly zero and yields a very presentable Pseudo Wigner-Ville Distribution. The 
ability to deemphasize these cross terms is a necessity if the desire is to evalute the energy 
change within a signal. For the remainder of this work, the window size selected (ie. N & M, 
N X M, or N by M) will be such that this criterion is met as closely as possible. 
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B. TIME AND FREQUENCY DEPENDANCE 

Since the criterion shown previously are time and frequency dependent, there must be 
a relation between the two that can be reduced to one parameter. The customary time and 
frequency relations of signal analysis are slightly changed in the context of the Wigner-Ville 
Distribution. Shown below are the time and frequency relations when working with the 
Wigner-Ville Distribution. 

U - DPxAt ; A f - - 2xDPxAf 

* * hnax 

where: t max = maximum signal length & f max = maximum frequency 

At = sampling time & Af = frequency resolution 
DP = number of data points 



The definitions of the size criterion discussed earlier for the smoothing window follow: 

_2 2 . 1 
® t^<i> ^ — 

4 

and since: 

o 2 t - (MAt) 2 - (NA(i>) 2 

and the time and frequency are related in the Wigner-Ville Distribution through the number 
of data points as shown: 



A(i> 



it 

2 DP At 
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Making this substitution into the size criterion, yields a very simple result for determining the 
size of the smoothing window required. 

DP 

MxN ± — 

71 

In looking at this result it seems as if the window is completely independent of time and 
frequency, but it must be rememberd that M and N are parameters that specify the lengths 
of the smoothing window in the time and frequency directions respectively. 

C. RELATION TO NUMBER OF DATA POINTS 

Another area of interest is the relation of the window size to the number of data points 
in the sampled signal. The best way to show this relation is through an example. For 
instance, a signal of a length 4.0 seconds that is sampled using 512 data points, will result in 
a At = 7.8125 msec and Af = 0.0625 Hz as calculated using the constraints of time and 
frequency as defined by the Wigner-Ville Distribution. Since the example has used 512 data 
points, the previous section has shown that chosing N = 10 and M = 16 will approximately 
satisfy the size criterion limitations. Additionaly, by definition of the smoothing window, the 
area of coverage will be: 

-2M<,m<L2M -2N<Ln<L2N 

This means that for the example shown, the coverage of the window will be 40 steps in the 
frequency domain and 64 steps in the time domain. 
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40xAf -2.5 



64 xA t - 0.5 



AREA -1.25 



where: Af = 0.0625 Hz & At = 7.8125 msec 

The AREA = 1.25 is a unitless quantity that defines the area that the smoothing window 
encompasses while smoothing one point in the Distribution. The unitless dimension is the 
result of the time being in seconds and the frequency in hertz. For comparison purposes, the 
total area of the distribution is the product of fm^and t max which is equal to 256.0, resulting 
in the window smoothing over 0.488% of the total Distribution, while smoothing one point. 

Now consider another example that uses 1024 data points. As shown earlier, the 
product of M and N must be greater than or equal to the number of data points divided by 
pi in order to meet the window size criterion as established by Reference 12. With respect to 
the earlier example, the size of the window (ie. N times M) must increase two times. 
Therefore, a selection of N = 10 and M = 32 will be used for this example. These values are 
very close to the required criterion and will suffice for this illustrative example. The signal 
length (t max ) will be 2.0 seconds resulting in f max = 256.25 Hz, At = 1.953 msec, and Af = 
0.125 Hz. In relation to the area of coverage of the window, the following results apply: 

40xA if -5.0 128xAt - 0.249984 AREA- 1.25 

where: Af = 0.125 Hz & At = 1.953 msec 

This is the same result as the 512 data point case, with the following exception. The total area 
of the Distribution (t max x f max = 512.0) has doubled, while the area of the smoothing window 
coverage has remained the same. This results in the window only covering 0.244% of the 
entire Distribution when smoothing a single point within the Distribution. In conclusion, the 
greater the number of sampled data points within the processed raw time signal, the smaller 
the area the smoothing window will encompass while smoothing the Distribution. An 
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important trade-off to consider, however, is the increase in computer time required to run 
the additional number of data points. 

D. TIME/FREQUENCY RESOLUTION 

Since the smoothing window used in this work gives the ability to vary the time and 
frequency size of the window, the following examples show the results that can be obtained. 
Figure 5, is the input signal that is sampled using 512 data points. It is a 400 Hz sine w'ave 
that has been computer generated with an amplitude of 2.0. In order to be consistent with the 
previous work, the size criterion established earlier will be used as a guide for the selection 
of M and N in the smoothing window. Figure 6, is the result of a window in which M and 
N are both equal to 13. Since this smoothing window does not require M and N to be equal, 
varying these parameters yielded some intersting results. Figure 7, shows the result of the 
input signal being smoothed with a 35 by 5 window. This means that N = 35 and M = 5, or 
in other words, the window is longer in the frequency domain than it is in the time domain. 
The resulting distribution has an amplitude that is much more constant throughout the entire 
time length, but the amplitude itself has been reduced considerably as compared to the 13 by 
13 example. Figure 8 is just a different view of this example, showing the loss of frequency 
resolution with a 35 by 5 window. Just the opposite case to this example is where the window 
used for smoothing is a 5 by 35 window. In this case, the window is longer in the time 
domain than the f requency domain and just the opposite results are obtained. Figure 9, shows 
the result of this window, while Figure 10 gives a view showing the degree of frequency 
resolution that can be obtained. This can be very beneficial in the monitoring of machinery, 
since in many cases the ability to differentiate between signals that are very close to one 
another is paramount. 
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400 HZ SINE WAVE 

SIGNAL LENGTH = 0.250 SEC / AMP = 2.0 
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Figure 5. 400 Hz Sine Wave Input 
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Figure 6. Pseudo Wigner-Ville Distribution (13 X 13 Window) 
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400 Hz SINE WAVE 
SMOOTHING WINDOW : 35 X 5 
AMPLITUDE = 2.0 




Figure 7. Pseudo Wigner-Ville Distribution (35 X 5 Window) 
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400 H» SINE WAVE 
FREQUENCY RESOLUTION : 35 X 5 
AMPLITUDE = 2.0 




r igure 8. Frequency Resolution (35 X 5 Window) 
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400 Hz SINE WAVE 
SMOOTHING WINDOW : 5 X 35 
AMPLITUDE = 2.0 




^S> ^cP 

< %y ^<51 



Figure 9. Pseudo Wigner-Ville Distribution (5 X 35 Window) 
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400 Hz SINE WAVE 
FREQUENCY RESOLUTION : 5 X 35 
AMPLITUDE = 2.0 
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igure 10. Frequency Resolution (5 X 35 Window) 
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E. TRANSITIONAL PEAK PHENOMENON 



In works that have given examples of the Wigner-Ville Distribution for transient 
signals, the point at which the signal’s frequency is changing from increasing to decreasing 
or vice versa has a large amplitude increase in the Wigner-Ville Distribution. For example, 
Rossano [Ref. 15] shows an artificially generated linear chirp and also an actual pump 
signature that has what appears to be a large concentration of energy at the transition point. 
In order to use the Wigner-Ville Distribution for transient machinery monitoring, there must 
be an explanation or way to quantify this "ghost peak" that is occuring at the transition point. 
To date, there has not been an explanation to this phenomenon, however, by looking at the 
unsmoothed Wigner-Ville Distribution, the reason for this is easily seen. Figure 1 1 is a linear 
chirp sine wave that is used as the input to the Wigner-Ville Distribution. Figure 12, shows 
the resulting unsmoothed Wigner-Ville Distribution of this linear chirp and as discussed and 
shown earlier, the cross-terms can be seen occuring exactly at the mid-point of the conflicting 
frequencies and times. These cross-terms build until they reach a maximum concentration 
at the transition or frequency change point. Since these are nothing more than cross-terms, 
to deemphasis these through the judicial use of the smoothing window. 
16 show the result of varying the smoothing window from a point of 
to entire elimination of the ghost peak. 



there must be a way 
Figures 13 through 
essentially no effect 
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UNEAR CHIRP SINE WAVE 
SIGNAL LENGTH = 0.256 SEC / AMP = 1.0 



Figure 11. 
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Linear Chirp Sine Wave Input 
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Figure 12. Unsmoothed Wigner-Ville Distribution 
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LINEAR CHIRP SINE WAVE 
SMOOTHED : 5 X 35 
AMPLITUDE = 1.0 







Figure 13. Pseudo Wigner-Ville Distribution (5 X 35 Window) 
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LINEAR CHIRP SINE WAVE 
SMOOTHED : 13 X 13 
AMPLITUDE = 1.0 







Figure 14. Pseudo Wigner-Ville Distribution (13 X 13 Window) 
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LINEAR CHIRP SINE WAVE 
SMOOTHED : 18 X 10 
AMPLITUDE = 1.0 
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LINEAR CHIRP SINE WAVE 
SMOOTHED : 35 X 5 
AMPLITUDE = 1.0 




Figure 16. Pseudo Wigner-Ville Distribution (35 X 5 Window) 
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This chapter has shown the benefits that may be gained through the use of the 
smoothing window. The frequency resolution and amplitude are the two major areas affected 
by the window changes. The last area presented, showed the reason for the large amplitude 
increase at the frequency transition point and how these "ghost peaks" may be smoothed and 
to what extent deemphasized by adjusting the smoothing window size. Within the next 
chapter, a thorough energy analysis of the discretized Pseudo Wigner-Ville Distribution will 
be conducted. 
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IV. WIGNER-VILLE DISTRIBUTION ENERGY ANALYSIS 



A. ENERGY DEFINITION 

As mentioned earlier, the Wigner-Ville Distribution is essentially the Fourier transform 
of the auto-correlation of an input signal. This may be thought of as the spectral density 
function in a three dimensional sense. The total energy of the Wigner-Ville Distribution was 
listed earlier as: 



where: 



1 

2n 




WDF{(S> , t) dcddt - t(r 



2 



i|/ 2 - mean square value of the signal 



There are a few points that must be made concerning this equation. The first is that the 
integration is from minus infinity to plus infinty and another is that it is the integration of 
the resulting Distribution of a real signal, not of an analytic signal as used in most discretized 
forms of the Wigner-Ville Distribution. Lastly is that this equation has been proven 
mathematically, which unfortunately does not take into account the fact that there are 
resulting cross-terms and that the Distribution is subsequently smoothed because of this. 
Since the last point is indicative of what is seen in actual real world problems, the energy must 
be defined in a discretized sense that accounts for the above mentioned points and also these 
unavoidable cross-terms. 
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B. EFFECT OF WINDOW SIZE ON ENERGY 



Since smoothing of the Distribution is a necessary requirement, the size of the applied 
window is something that must be decided upon. Another important point is knowing what 
effect this will have on the energy of the Distribution. In the previous chapter, the size 
criterion that was established by Cartwright [Ref. 19] and later used by Wahl and Bolton in 
Reference 12 was discussed. Within this discussion, the visual benefits of varying the window 
size in the time and frequency domain were presented. In order to define the energy, the 
effect on the energy of the Pseudo Wigner-Ville Distribution as caused by varying these 
windows will be presented. 

Figure 17 is a 400 Hz sine wave of length 0.256 seconds and amplitude 2.0 that will be 
sampled and have the Wigner-Ville Distribution calculated. Figure 18 is the resulting 
smoothed Wigner-Ville Distribution that has been processed using a 13 by 13 sized smoothing 
window. Figure 18 also shows that the volume is equal to 0.0709118 V 2 . This is 
representative of the energy contained within the Wigner-Ville Distribution of the signal in 
Figure 17 and is nothing more than the volume under the resulting Distribution. By 
definition, this volume should be equal to the mean square value of the 400 Hz sine wave, 
however, due to the points mentioned above, this is not going to hold true and will be 
discussed in depth later. Figures 19 and 20 are also smoothed Distributions of this same 400 
Hz signal only using a 35 by 5 window in the case of Figure 19, and a 5 by 35 window in 
Figure 20. Since N times M (the size criterion) are very close in these three examples, the 
resulting energy values of the three cases are essentially the same. The point to be emphasized 
from this is that as long as the area of the window does not vary significantly, the energy will 
be represented in a constant manner and independent of the time and frequency window size 
selected. The last example that will be presented in this section is that this result is 
independent of the frequency of the signal as shown by the combined 100 and 400 Hz sine 
wave example of Figure 2 1 , in which the energy is equal to 0. 1 4 1 83 1 4 V 2 . This is essentially 
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twice the value of the single sine waves of Figures 18 through 20. Of special importance for 
transient machinery applications is Figure 22, which is a transient or ramped sine wave that 
shows the same resulting energy as the previous stationary signals. This shows that the 
Wigner-Ville Distribution energy is independent of the nature of the signal allowing for the 
possibility of monitoring a piece of machinery during both it’s steady state and transient 
phases of operation. 
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400 HZ SINE WAVE 

SIGNAL LENGTH = 0.256 SEC / AMP = 2.0 




O'Z 0't 00 01- 0 3- 

(sjioa) gamndiAiv 



Figure 17. 400 Hz Sine Wave Input 
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Figure 18. Pseudo Wigner-Ville Distribution Energy (13 X 13 Window) 
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400 Hz SINE WAVE 
SMOOTHING WINDOW ~as x 5 
AMPLITUDE = 2.0 




Figure 19. Pseudo Wigner-Ville Distribution Energy (35 X 5 Window) 
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400 Hz SINE WAVE 
SMOOTHING WINDOW ; 5 X 35 
AMPLITUDE = 2.0 




Figure 20. Pseudo Wigner-Ville Distribution Energy (5 X 35 Window) 
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Figure 21. Combined Sine Wave Distribution Energy (13 X 13 Window) 
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RAMPED SINE WAVE 
SMOOTHING WINDOW : 13 X 13 
AMPLITUDE = 2.0 
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Figure 22. Ramped Sine Wave Distribution Energy (13 X 13 Window) 
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The previous figures show that if the Wigner-Ville Distribution is going to be used to 
accurately portray or maintain signal information, then the window size must at least remain 
constant and better yet to keep the identical window size for each measurement of 
comparison. This requirement to maintain a constant sized smoothing window is not the only 
constraint that must be meet for energy analysis purposes. The next section will discuss 
another of these elements. 

C. TIME DEPENDENCE OF WIGNER-VILLE ENERGY 

Not only is the window size going to influence the energy result, but the time length 
of the signal will also play a large part. This varies considerably from the conventional 
knowledge of the mean square value of a signal being independent of time. The standard 
mean square value will be the same for a signal of time length 2.0 seconds as it will for a 
signal of 4.0 seconds. For instance, the mean square value of a sine wave x(t) is: 



where: x(t) = X sin (cot) 

To compute the mean square value of this sine wave, it is independent of time and will equal 
0.5 for X = 1.0. This is not the case for the mean square value or total energy of the signal 
as found using the discretized form of the Wigner-Ville Distribution. The following figures 
will show this second important point that must be realized when comparing Wigner-Ville 
energies of various signals. Figure 23 is the 400 Hz input signal of length 0.128 seconds. This 
time length is half of that of the 400 Hz signal used in the previous section. The resulting 
Wigner-Ville Distribution is shown in Figure 24 with an energy value of 0.0354558 V 2 . As 
can easily be seen, this resulting energy value is not the same and in fact, the energy of Figure 




o 
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24, is one half that of the earlier 400 Hz signals of length 0.256 seconds. Had the signal been 
of a length 0.512 seconds or twice as long as the earlier examples, the resulting Wigner-Ville 
energy would have doubled. The window size used for this Pseudo Wigner-Ville Distribution 
was a 13 by 13 sized window. 



44 



400 HZ SINE WAVE 

SIGNAL LENGTH = 0.128 SEC / AMP = 2.0 




(siioa) ganindiAiv 



Figure 23. 400 Hz Sine Wave (T^ — 0.128 sec) 
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Figure 24. Pseudo Wigner-Ville Distribution Energy (T L = 0.128 sec) 
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This has yielded a second concern that must be kept in mind when using the energy of 
the discretized Pseudo Wigner-Ville Distribution. Since the energy is proportional to the 
signal length, an easy conversion can be made by using the lengths of the individual raw input 
signals. This time dependence holds true for not only steady state signals as shown, but also 
for signals of a transient nature such as previously shown with Figure 22. Throughout these 
last two sections, the comparison of energy values of signals has been mentioned a great deal. 
The thought process behind this is that if there is an amplitude change in the input signal, this 
will be reflected in the resulting Pseudo Wigner-Ville Distribution. By ensuring that these 
previous two conditions are met, an investigation into the resulting energy changes from an 
amplitude increase or decrease of an input signal can be made. 

D. ENERGY PROPORTIONALITY OF THE DISTRIBUTION 

The mean square value of a simple sine wave is the amplitude squared divided by two. 
This was shown in an earlier equation and now will be used to show the proportionality of the 
Wigner-Ville Distribution. For instance, the following applies for signals of differing 
amplitudes as shown: 

amplitude - 1.0 i|r 2 - 0 . 5 

amplitude - 2.0 i^/ 2 — 2.0 

This amplitude increase from 1.0 to 2.0 represents an increase of four in the mean square 
value of the signal. Figures 25 and 26 are the raw input sine waves of amplitude 1.0 and 2.0 
respectively. The resulting Pseudo Wigner-Ville Distributions are shown in Figures 27 and 
28, which acurately portray the four-fold energy increase. These Distributions have been 
smoothed using a window of 13 by 13. 
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400 HZ SINE WAVE 

SIGNAL LENGTH = 0.256 SEC / AMP = 1.0 




Figure 25. 400 Hz Sine Wave (Amplitude = 1.0) 
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Figure 26. 400 Hz Sine Wave (Amplitude = 2.0) 
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Figure 27. Pseudo Wigner-Ville Distribution (Amplitude = 1.0) 
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400 Hz SINE WAVE 
SMOOTHING WINDOW ; 13 X 13 
AMPLITUDE = 2.0 
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Figure 28. Pseudo Wigner-Ville Distribution (Amplitude = 2.0) 
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These previous four figures show the proportionality of the energy within the Wigner- 
Ville Distribution. These two cases are just an example to show the proportionality of the 
energy. Runs were made for many more cases (ie. X = 4.0, X = 8.0, etc) and the 
proportionality relation held true as shown with the previous Figures. This was another 
required element in order to use the energy as the basis for a method of machinery 
monitoring. Had the cross-terms increased in such a fashion as to not allow the Distribution 
to show the proportional energy change, the energy content of the Pseudo Wigner-Ville 
Distribution would have been rendered useless for the purposes of comparison and evaluation. 
The next chapter will use these requirements and basics that have been developed to show an 
actual application of how the Wigner-Ville energy may be used to conduct machinery 
monitoring. 
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V. GEAR MODEL APPLICATION 



A. DESCRIPTION OF GEAR MODEL 

The simple complexity gear model used in this work was built based upon a model used 
by D.K. Carlson [Ref. 20] for Artificial Neural Network applications for machinery 
monitoring. It may be of a simple nature, however, it provides all of the necessary 
components that may be found in more complex systems. A schematic of the gear model can 
be seen in Figure 29. The two major changes of this model over the prototype used by 
Carlson is that the components are much larger enabling a better signal to noise ratio and 
secondly, the motor used to drive the gear model is mounted on a separate base to reduce 
motor enduced vibrations. 

The model consists of a single reduction gear train, with the pinion gear being a 90 
tooth spur gear and the driven gear a 120 tooth spur gear. The gears were manufactured by 
Martin Corporation with a 1/2 inch bore and a 14.5 degree pressure angle. Aluminum blocks 
housed the NICE 1/2 inch bore radial ball bearings that were used to support the shaft and 
gears. These aluminum blocks were bolted to the support base which was a 1.0 inch thick 
plexiglass slab. Driving the single shaft is a General Electric 1 /6 th HP motor, that is 
controlled by the use of a Bodine Electric Company combination rectifier and variable 
potentiometer speed controller. The shaft speed was monitored through use of a Power 
Instruments Model 1720 RPM Optical Proximeter. [Ref 20] 

B. SIGNAL SAMPLING EQUIPMENT 

This section will describe the components shown in Figure 30, that were used to sample 
the vibrational signal of the gear model. The signal was monitored using an ENDEVCO 
model 303A03 accelerometer that was amplified and powered with a PCB model 483B07 
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power unit. The signal was then processed through a Krohn-Hite model 3342 analog filter 
configured in a band pass mode, prior to being sent to the HP 3565S system for sampling and 
storage. Lastly, the signal was transfered to the VAX 3520 workstation for the calculation and 
plotting of the Pseudo Wigner-Ville Distribution. 

1. ENDEVCO Model 303 A03 Accelerometer 

The ENDEVCO Isotron PE accelerometer is a minature accelerometer that has a 
medium range high frequency capability. The operation of the accelerometer is based upon 
a piezoelectric quartz transducer sensing element. The sensitivity of the accelerometer is 10 
mV/g with a resonant frequency of 70 kHz. The resolution available is 0.02 g, with a 
maximum range of + 500 g. 

The accelerometer was mounted on top of the aluminum bearing housing that 
supports the shaft and specifically at the position shown in Figure 30. The accelerometer wass 
attached by means of a mounting wax which yields a maximum operating range out to 15k - 
20k Hz. This accelerometer signal was then sent through the power unit described below. 

2. PCB Model 483307 Power Unit 

This power unit is used to power the low impedance quartz transducers and 
amplify the signal if desired. The power unit provides an adjustable 2 to 20 mA current for 
purposes of transducer excitation. The gain adjustment is available from 0 to 100 and is set 
through the use of a ten turn vernier gain pot and a three position gain multiplier switch. For 
the purposes of these tests, the gain was set to 20 before sending the signal to the analog filter. 

3. Krohn-Hite Model 3342 Anolog Filter 

The model 3342 variable filter is a digitally tunned filter that will function as a 
High-Pass or Low-Pass filter. When the two channels of the filter are connected together, the 
filter will function as a band pass filter, which is the configuration for the gear model work. 
The range of the filter is from 0.001 Hz to 99.9 kHz as adjusted by three rotary decade 
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switches and a rotary six positiohn multiplier switch. In addition, the filter unit has a gain 
setting of unity (0 dB) and 10 (20 dB). The signal was not further amplified using this 
capability of the filter prior to sampkling by the HP 3565S system. 

4. HP 3565S Measurement Hardware System 

The HP 3565S measurement hardware with the HP- Vista sof tware uses a HP-9000 
series computer system for controlling purposes. The HP 3565S system is a modular 
multichannel system that can process data in both the time and frequency domain. The 
system is capable of handling 64 source and input modules, but is presently configured for 
eight. The HP-Vista software allowed for the sampling and storage of the preprocessed 
accelerometer signal that was later processed through the Pseudo Wigner-Ville 
distribution. 
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Figure 29. Gear Model Schematic 




Accelerometer 



Figure 30. Signal Sampling Equipment Schematic 
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C. ESTABLISHING MONITORING FREQUENCIES 



This section will describe how the monitoring frequencies have been determined and 
how they will be used to monitor the gear model and in future applications a piece of 
transient machinery. As shown earlier in Chapter III, the energy contained in a signal is 
accurately portrayed within the resulting time-frequency domain of the Wigner-Ville 
Distribution. Through the Wigner-Ville Distribution energy that is contained in the 
established monitoring frequencies, faults may be detected. The shaft frequency is just the 
speed of the shaft, while the gear mesh frequency is the number of teeth on the rotating gear 
multiplied by the shaft speed or frequency. Both of these frequencies have associated 
harmonics that can be of equal importance as the frequencies themselves in the monitoring 
of machinery, It is for these reasons that the frequencies listed below have been chosen based 
upon a shaft frequency of 15 Hz. 

1. Monitoring Frequencies 

For the purposes of this work, just the shaft and gear mesh frequencies of the 
model will be included. Just as in standard established machinery monitoring methods, the 
shaft frequency and associated harmonics can provide a great deal of information for fault 
detection and diagnosis. Along with the shaft frequencies, the gear meshing frequencies 
provide valuable information of faults and will also be monitored. The purpose of the two 
examples will be to show that the Wigner-Ville Distribution energy (ie. volume) changes as 
a result of damage and will prove to be a very valuable tool in future machinery monitoring 
methods. The established monitoring frequencies have been filtered as follows: 

• Shaft Input 1: (5 - 100 Hz) This is the spectrum of signals in the lower frequency range 
selected to determine the total energy contained within the shaft frequency and it’s 
harmonics. 

• Shaft Input 2: (10 - 20 Hz) This frequency range captures just the shaft frequency. 

• Shaft Input 3: (25 - 35 Hz) This frequency range captures the 1 st shaft harmonic. 
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• Shaft Input 4: (40 - 50 Hz) This frequency range captures the 2 nd shaft harmonic. 

• Shaft Input 5: (55 - 100 Hz) This frequency range captures the upper harmonics of the 
shaft frequency. 

• Gear Mesh Input 1: (1250 - 1450 Hz) This frequency range captures the gear meshing 
frequency. 

• Gear Mesh Input 2: (2600 - 2800 Hz) This frequency range captures the 1 st gear mesh 
harmonic. 

• Gear Mesh Input 3: (3950 - 4150 Hz) This frequency range captures the 2 nd gear mesh 
harmonic. 

2. Combined Monitoring Technique 

Using the Wigner-Ville energy as an input to an Artificial Neural Network is a 
follow-on goal of this research. It is for this reason, that the above established frequencies 
have been identified as inputs. Within this scheme, the Wigner-Ville is used to characterize 
the signal and give a single energy value at the particular frequency of interest as an input to 
an Artificial Neural Network. By filtering these signals prior to processing, just the 
frequency of interest will be captured along with it’s respective Wigner-Ville energy content. 
The Wigner-Ville energy data can now be used to produce a comprehensive training set for 
the Artificial Neural Network. By pre-processing the inputs to a Neural Network with the 
filtering and the Wigner-Ville Distribution, the number of inputs is greatly reduced. Since 
these associated energy levels are fairly stable for a given machinery condition, the Neural 
Network will be capable of detecting patterns in energy level shifts that are associated with 
particular machinery faults. To support this last statement, the following section will show 
two examples of machinery faults imposed on the gear model and the response of the Wigner- 
Ville energy to these faults. 



58 



D. MACHINERY FAULT EXAMPLES 



1. Mechanical Looseness 

This fault was imposed on the gear model by loosening the aluminum support 
blocks that house the bearings and support the shaft. Specifically, the blocks that support the 
driver shaf t were loosened. The degree of looseness imposed was difficult to judge and made 
this particular fault a difficult one to diagnose. The blocks were not completely undone, but 
loosened just a small amount from their original firmly bolted posture. Shown in the table 
below are the results of averaging two good runs and two post-damage runs and comparing 
the respective energy values at the different frequencies. These energy values are nothing 
more than the volume under the Pseudo Wigner-Ville Distribution as discussed earlier. Only 
the shaft frequencies are shown for purposes of clarity and since it is in these monitored 
frequencies that a fault of mechanical looseness should present itself. 
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Table I. MECHANICAL LOOSENESS 



MECHANICAL LOOSENESS 

ENERGY VALUES (V~2) 



Frequency 
Shaft Spectrum 
Shaft Frequency 
1st Harmonic 
2nd Harmonic 
Upper Spectrum 



Good 


Damaaed 


41.0 


28.0 


1.07 


0.83 


2.31 


1.22 


1.23 


2.40 


31.4 


26.31 



As can be seen in the table above, the energy content approximately doubled at the second 
shaft harmonic. This is a very encouraging result for the purposes of using the Pseudo 
Wigner-Ville Distribution energy as a fault detection tool. 

2. Gear Tooth Damage 

This particular fault was imposed on the gear model by shaving a gear tooth on 
the driver gear. The degree of damage was much easier to predict in this case and made this 
particular fault an easier one to control. The gear tooth that was damaged had approximately, 
1 /8 th of the leading edge removed. Shown in the table below are again the results of 
averaging two good runs and two post-damage runs and comparing the respective energy 
values at the different frequencies. For this case, just the frequencies associated with the gear 
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mesh are presented, since it is in these monitored frequencies that a fault involving the 
damage of a gear train would present itself. 

Table II. GEAR TOOTH DAMAGE 



GEAR TOOTH DAMAGE 



ENERGY VALUES 


C\J 

< 

> 


MESH FREQUENCY 


GOOD 


DAMAGED 


Fundamental 


17.14 


154.4 


1st Harmonic 


0.53 


3.44 


2nd Harmonic 


0.16 


0.48 



As shown in the above table, the energy content greatly increased at the fundamental gear 
mesh frequency. Additionaly, there was also increases at the harmonics and shown in Figures 
31 and 32 is a graphical presentation of the energy increase that occured due to the damaged 
tooth. 



61 



GEAR TOOTH DAMAGE 
GEAR MESH INPUT ONE 



FILTERED 1250 - 1450 HZ 




Figure 31. Gear Mesh Frequency (Pre-Damage) 
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GEAR TOOTH DAMAGE 
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Figure 32. Gear Mesh Frequency (Post-Damage) 
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This chapter has shown two examples of machinery damage and the resulting changes 
of the Pseudo Wigner-Ville Distribution energy. By applying the window and energy 
requirements discussed in Chapters III and IV, the Wigner-Ville energy has yielded a definite 
change to allow a means for machinery monitoring. The Wigner-Ville energy will make a 
stable input to a Neural Network and provide a quick efficient method by pre-processing the 
machinery signal with the Wigner-Ville Distribution as shown. This pre-processing will allow 
the number of inputs to the Neural Network to remain at a reasonable level and provide for 
a manageable, quicker hardware system. 
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VI. CONCLUSIONS 



The Pseudo Wigner-Ville Distribution has been shown to provide an excellent 
characterization of an input signal. Based upon the research conducted, it will also provide 
a means for machinery condition monitoring and diagnostics. The research has also shown 
that in order to use the Wigner-Ville Distribution in a machinery monitoring scheme, the 
following conditions or conclusions must be met and understood. 

• The smoothing window size applied must remain constant for purposes of comparing 
various signals. 

• The individual time and frequency lengths of the smoothing window may be adjusted 
to provide a better frequency resolution or constant amplitude. 

• The energy contained within the Pseudo Wigner-Ville Distribution of a signal is 
proportionaly time dependent. 

• The energy of the Distribution is proportional to the mean square value of the input 
signal and is not adversely affected by the cross-terms. 

• Lastly, the damaged gear model conditions showed a definite change in the energy of 
the Pseudo Wigner-Ville Distribution at the respective monitoring frequencies. 
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VII. RECOMMENDATIONS FOR FUTURE RESEARCH 



The Pseudo Wigner-Ville Distribution computer program that was established here at 
the Naval Postgraduate School is very stable. The pronounced peaks at the transition points 
are now definable and additionally, with the energy of the Pseudo Wigner-Ville Distribution 
now being a definable and understood quantity, the next additional research in this area may 
include: 

• Establish a baseline or data bank of "good" Wigner-Ville Distribution data for the 
present gear model. 

• Conduct additional runs of various levels of damage of gear model components for 
comparison with good data. 

• Establish a Torpedo Ejection Pump baseline of good data as with the gear model. 

• Continue with the research of linking the Pseudo Wigner-Ville Distribution as an input 
to a Neural Network Application. 
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APPENDIX COMPUTER PROGRAMS 



The following programs have all been developed here at the Naval Postgraduate School 
and used in the conduct of the research presented. The language used is FORTRAN 77 with 
the plotting programs using calls to CA-DISSPLA software. 



A. WVD3 COMPUTER PROGRAM 



PROGRAM WVD3 

* 

************************************************************ 



* THE FOLLOWING INDIVIDUALS HAVE BEEN INSTRUMENTAL IN THE * 

* DEVELOPMENT OF THIS PROGRAM AND PREVIOUS VERSIONS: * 



★ ★ 

★ * 

* Graham Rossano * 

* * 

* Jose G. Calahorrano * 

* * 

* Scott G. Spooner * 

* * 

* Prof. Y.S. Shin * 

★ ★ 

* * 



************************************************************ 

* 

* This program calculates the Wigner Distribution 

* Function of a signal. It includes the option of applying a 

* smoothing function and reducing the distribution to different 

* sizes while writing the results to files for subsequent 

* plotting. 

* 



************************************************************ 



VERSION 3 



************************************************************ 



VARIABLE DECLARATION 



****★***★*★****★★*★★★*****★*★★*★★*****★*****★*★★*★***★★***★* 



★**★*★★★★★... setting number of data points — ★***★★★★***★★ 



integer dp,dp2,mvopt ,redopt ,dimt,dimf ,mm,nn,rf ,rdp,rdp2, 

: nf ,mt,nf2,mt2 

PARAMETER (DP = 512) 

real pi , t in(dp) ,ain(dp) , rawa(dp) , rwdf (dp*2,dp) , rswdf (dp*2,dp) , 
: rawt(dp),dtSLin,delt,dt,asLin,meanv,mtime,del1 ,del2, 

: const , t , sum, sumb, val , sval , wdf (dp*2,dp) ,coef ,df ,AREA 

complex s(dp*2),dun,c(dp*5),duiTL5,ckjn2 
character*25 inname 

c ommon/ i n/ mvopt , r edopt , d i mt , d i mf , mm, nn , r f , r dp , rdp2 , 
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nf ,mt,nf2,mt2 

coomon/rl/ pi ,t in,ain,rawa, rwdf ,rswdf , 

: rawt,dtsun,delt ,dt,asLm,meanv,mtime,del1 ,del2, 

: const, t ,sun,sunb,val ,sval ,wdf ,coef ,df ,area 

common/comp/ s,duri,c,dun3,dLin2 
dp2 = dp* 2 

* 

************************************************************ 



* INITIALIZING VARIABLES * 

************************************************************ 



— Print description of program — 
print* 

print*, 1 Program to calculate the' 
print*,' Pseudo Wigner-Ville Distribution' 
print* 



— Set initial values — 



print*,' Enter name of signal input file' 

read(5,904) inname 

print* 

print*,' Do you wish to remove the mean value?' 
print*,' Enter 1 for yes - 0 for no' 
read(5,910) mvopt 
print* 

print*,' Input the reduction size desired' 
print*,' input 1 for 64 by 32 ' 

print*,' input 2 for 128 by 64 ' 

print*,' input 3 for 128 by 128 ' 

print*,' input 4 for 256 by 128 ' 

read(5,910) redopt 
print* 



print*,' Input the smoothing window size desired' 
print*,' input the size for the frequency axis' 
read(5,9l0) nf 

print*,' input the size for the time axis' 
read<5,910) mt 

* 

* 

* 

************************************************************ 
* * 

* calculation of some common used constants * 

* * 

************************************************************ 



pi = 4.0 * atan(I.O) 



**************************************************************** 

* THIS IS THE CALCULATION PART OF THE PROGRAM * 

**************************************************************** 

open(4,f i le=inname,status='old' ) 
rewind 4 

* PRELIMINARY CALCULATIONS: 

CALL DATAIN3 

write(6,*> , ' done with datain3' 

CALL DTCALC3 
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write(6,*) , 1 done with dtcalc3' 

CALL MEAN3 

write(6,*) , ' done with mean3' 
write(6,*) , ' mean value removed = 1 , meanv 

* SIGNAL MODIFICATIONS: 

* Window application (modified hamming window) 

CALL HAMMG3 

write(6,*) , ' done with hafrmg3' 

* Conversion of real signal to an analytical one 

CALL ANAL3 

write(6,*) , 1 done with anal3' 



* CALCULATION OF THE WIGNER VILLE DISTRIBUTION: 
CALL WIGH3 

write<6,*) , 1 done with wigh3‘ 



* REDUCTION & SMOOTHING OF THE WIGNER VILLE DISTRIBUTION: 

CALL REDUCE3 

write(6,*) , ' done with reduce3' 

CALL SMOOT H3 

write(6,*) , ' done with smooth3‘ 

*********************************************************** 

* * 

* WRITING OF THE DIFFERENT ARRAYS GENERATED * 

* TO FILES FOR PLOTTING PURPOSES. * 

* * 

*********************************************************** 

* 

* OPENING OF THE OUTPUT FILES 

* 

OPEN (UN I T=7, F I LE= ' RAWT I ME .OUT ' , STATUS= ' NEW 1 ) 
OPEN(UNI T=8, F ILE= ' WNDWT IM.OUT ' , STATUS= 1 NEW' ) 

OPEN (UN I T=9 , F I LE= ' RWDF . OUT ■ , STATUS= ' NEW ' ) 
OPEN(UN1T=10, F I LE= 'RSWDF .OUT 1 > STATUS =, NEW' ) 
0PEN(UNIT=1 1 , FILE='WDF.OUT ' , STATUS 1 'NEW' ) 

• 

* WRITING OF RAW AND WINDOWED TIME SIGNAL TO FILES 

* 

DO 300 I = 1,DP 

WRITE(7,930) TIN( I ),RAWA( I ) 

WRITE(8,930) TIN( I ),AIN( I ) 

300 CONTINUE 
* 

* WRITING OF REDUCED/UNSMOOTHED WVD TO FILE 

* 

DO 400 I = 1, DP/mm 
DO 400 J = 1 ,DP2/nn 
WRITE(9,906) RWDF(J,I) 

400 CONTINUE 

* 

* WRITING OF REDUCED & SMOOTHED WVD TO FILE 
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00 500 I = 1 ,DP/nri 
DO 500 J = 1 ,DP2/nn 

WRITE( 10,906) RSUDF(J.I) 

500 CONTINUE 

* 

* WRITING OF ENTIRE UDF TO A FILE 

* 

DO 600 I = 1,DP 

DO 600 J = 1.DP2 

URITE(11,906) UDF(J,I) 
600 CONTINUE 

* Format statements 

904 format(a25) 

906 format(2X,g16.8) 

910 formate i6) 

930 format(2x,g16.8,5x,g16.8) 

END 

* SUBROUTINES 

include 'ANAL3. INC' 
include 'CORR3.INC' 
include 'DATAIN3. INC 1 
include 'DTCALC3. INC 1 
include 1 FFT3. INC 1 
include 1 HAMMG3 .INC 1 
include 'HEAN3.INC 1 
include 'SHOOTH3. INC 1 
include 'REDUCE3. INC' 
include 'WIGH3. INC' 



B. SUBROUTINES OF WVD3 



Listed below in alphabetic order are the subroutines used in the main program WVD3. 



SUBROUTINE ANAL3 

* 

********************************************************** 



* * 

* This subroutine converts a real signal to an * 

* analytic one * 

* * 



********************************************************** 



* 

integer dp,dp2,mvopt , redopt,dimt # dimf f rmi # nn, rf , rdp f rdp2, 

: nf ,mt,nf2,mt2 

PARAMETER (DP = 512) 

real pi # tin(dp),ain(dp),rawa(dp),rwdf (dp*2,dp), rswdf (dp*2,dp), 
: rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1 ,del2, 

: const,t,sum,sumb,val ,sval ,wdf (dp*2,dp),coef ,df ,area 

complex s (dp* 2), dun, c(dp*5) # duin3 , dc»n2 
character*25 inname 

* 

common/ i n/ mvopt , redopt , di mt , di mf , rmi # nn, rf , rdp, rdp2 f 
: nf ,mt,nf2,mt2 

common/ rl/ pi ,tin,ain,rawa,rwdf , rswdf , 

: rawt,dtsum,delt,dt,asum,meanv,mtime # del1 ,del2, 

: const , t , sum, sumb, va l , sva l , wdf , coef , df , a rea 

common/comp/ s,dum,c,dum3,ckjn2 
dp2 = dp* 2 

* 

* 

do 100 i=1,dp 
sum=0.0 

do 200 j = 1 , dp 
sumb=0.0 

if(i-j.eq.O) go to 200 



n=i-j 

val=pi*n/2. 

sval=sin(val ) 

sumb=ain( j )*sval*sval/val 

200 sum=simsunb 

s(i )=cmplx(ain( i ),sum) 

100 continue 

END 
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SUBROUTINE DATAIN3 

************************************************************** 

* * 

* This subroutine assures the input file is in the following * 

* format: time amplitude (2x,e16.8,5x,e16.8) with an end of * 

* file indicator of a last entry 9999., 9999. * 

* * 
**#**#**t********t*ttttttttttttt**tttttttttttttt*ttt#t*<r***t*t 
* 

★ 

★ 

integer dp,dp2,mvopt ,redopt,dimt,dimf ,mm,nn, rf ,rdp,rdp2, 

: nf ,mt,nf2,mt2 

PARAMETER (DP = 512) 

real pi , t in(dp),ain(dp), rawa(dp), rwdf (dp*2,dp), rswdf (dp*2,dp), 
: rawt(dp),dtsum,delt ,dt,asum,meanv,mtime,del1 ,del2, 

: const , t , sum, sumb, va l , sval , wdf (dp*2,dp) , coef ,df ,AREA 

complex s(dp*2),duTi,c(dp*5),dum3,dum2 
character*25 inname 

* 

common/ i n/ mvopt , redopt ,dimt ,dimf ,mm, nn, rf , rdp, rdp2, 

: nf ,mt,nf2,mt2 

common/rl/ pi , tin, a in, rawa, rwdf , rswdf , 

: rawt,dtsun,delt ,dt,asun,meanv,mt ime,del 1 ,del2, 

: const, t,sun, sumb, va l, sval, wdf , coef ,df, area 

common/comp/ s,dun,c,dun3,dum2 
dp2 = dp* 2 

★ 

********simple loop to read in time & amplitude************* 
do 100 j=1 ,dp 

read(4,902) rawt(j) , rawa(j) 
tin(j) = rawt(j) 
ain(j) = rawa(j) 





100 


continue 




★ 


FORMAT STATEMENT 




902 


FORMAT (2X,E16.8,5X,E16.8) 




200 


END 
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SUBROUTINE DTCALC3 



********************************************************** 

* * 

* This subroutine calculates the delta t of the signal * 

* * 

********************************************************** 

* 

* 

* 

integer dp,dp2,mvopt , redopt ,dimt,dimf ,mm,nn,rf , rdp,rdp2, 

: nf ,mt,nf2,mt2 

PARAMETER (DP = 512) 

real pi ,tin(dp),ain(dp),rawa(dp),rwdf(dp*2,dp),rswdf(dp*2,dp), 
: rawt( dp), dtsum, del t ,dt ,asum,meanv,mtime,del1 ,del2, 

: const, t,sum,sumb,val ,sval ,wdf (dp*2,dp),coef ,df ,AREA 

complex s(dp^2),dtin,c(dp*5),dum3,dum2 
character*25 inname 

* 

common/ i n/ mvopt , redopt , di mt , d i mf , mm, rm, rf , rdp, rdp2 , 

: nf ,mt,nf2,mt2 

common/rl/ pi , t in,ain, rawa, rwdf , rswdf , 

: rawt,dtsum,delt ,dt ,asum,meanv,mt ime,del 1 ,del2, 

: const, t,sum,sumb,val ,sval,wdf ,coef ,df ,area 

common/comp/ s,dum,c,dum3,dum2 
dp2 = dp* 2 

* 

* 

************************************************************ 

dtsum = 0.0 

do 100 i = 1 , dp-1 

delt = t i n( i+1 ) - tin(i) 
dtsum = dtsum + delt 

1 00 cont i nue 

dt = dtsum / (dp-1 . ) 



END 
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SUBROUTINE FFT3 



********************************************************* 

* 

This subroutine calculates the Fast Fourier Transform * 

* 

********************************************************* 



integer dp,dp2,mvopt, redopt,dimt,dimf ,mm,nn, rf ,rdp,rdp2, 

: nf ,mt,nf2,mt2 

PARAMETER (DP = 512) 

real pi , t i n(dp) , ai n(dp) , rawa(dp) , rwdf (dp*2,dp) , rswdf (dp*2,dp) , 
: rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1 ,del2, 

: const, t ,sum,sumb,val ,sval ,wdf (dp*2,dp),coef ,df 

complex s(dp*2),dum,c(dp*5),durt3,ckin2 
character*25 inname 

common/ in/ mvopt,redopt,dimt,dimf ,mm,nn,rf ,rdp,rdp2, 

: nf ,mt,nf2,mt2 

common/rl/ pi , tin, ain,rawa, rwdf, rswdf , 

: rawt,dtsum,del t,dt,asum # meanv,mtime,del1 ,del2, 

: const, t ,sum,sumb,val ,sval ,wdf ,coef ,df ,area 

common/comp/ s,dum,c,cki , n3,dun2 
dp2 = dp*2.0 



const=dp2 

val=alog(const )/alog(2. )♦. 1 

j = 1 

do 40 i=1 ,dp2-1 

if (i .ge. j ) go to 10 
dun3=c( j ) 
c( j )-c( i ) 
c ( i )=dun3 
10 k=dp 

20 if (k.ge.j) go to 30 

j=j-k 
k=k/2 
go to 20 
30 j=j+k 

40 continue 

do 70 n=1,val 
coef=2**n 
coef=coef/2 
dum2=cmplx(1 . ,0. ) 

ckin=cmp l x( cos ( p i /coef ) , - s i n( p i /coef ) ) 

do 60 j=1 ,coef 

do 50 i=j,dp2,2*coef 
i i=i+coef 
dum3=c(i i )*dum2 
c(i i)=c(i)-dum3 
c(i )=c( i )+dun3 
50 continue 

dum2=dum2*dum 
60 continue 
70 cont i nue 

END 
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SUBROUTINE HAMMG3 



******************************************************** 

* * 

* This subroutine applies a modified hamming window * 

* to the signal ain(t) * 

* ★ 
******************************************************** 

* 

* 

integer dfc>,dp2,mvopt,redopt,dimt,dimf ,mm,m,rf ,rdp,rdp2, 

: nf ,mt,nf2,mt2 

PARAMETER (DP = 512) 

real pi , tin(dp),ain(dp), rawa(dp) , rwdf (dp*2,dp), rswdf (dp*2,dp), 
: rawt(dp) # dtsun # delt,dt # asun,rr>eanv,mtifne,del1 ,del2, 

: const,t,sun,sumb,val,sval ,wdf <dp*2,dp) ,coef ,df ,area 

complex s(dp*2),ckin,c<d£*5),dun3,dcin2 
character*25 imame 

* 

common/ in/ rnvopt , redopt ,dimt ,dimf ,ntn,nn,rf , rdp,rdp2, 

: nf ,mt # nf 2,mt2 

common/rl/ pi , tin, ain,rawa, rwdf , rswdf , 

: rawt,dtsun,delt,dt,asun,meanv,mtime,del1 ,del2, 

: const, t ,sun # sumb,val ,sval ,wdf ,coef ,df ,area 

common/comp/ s,dun,c ,dum3,dum2 
dp2 = dp* 2 

* 

mt ime=dp*dt 
del1=0.1*mtime 
del2=0.9*mt ime 
cons t=pi /dell 

do 100 j = 1 , dp 

t = (j-1) * dt 

if (t. le.dell) then 

ain(j) = ain(j) * ( .5*0 . -cos(const*t))) 
elseif ((t ,ge.del2) .and. (t . le.mtime)) then 
ain( j)=ain( j)*(.5*(1 . -cos(const*(mtime-t )))) 
else 
end if 

100 continue 
END 
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SUBROUTINE MEAN3 



* 

********************************************************** 
* * 

* This subroutine calculates and removes the mean value * 

* of the signal . * 

* * 
********************************************************** 



integer dp,dp2,mvopt , redopt,dimt,dimf ,mm,nn, rf , rdp, rdp2, 

: nf # mt , nf2,mt2 

PARAMETER (DP = 512) 

real pi ,tin(dp),ain(dp),rawa(dp),rwdf (dp*2,dp), rswdf (dp*2,dp), 
: rawt(dp),dtsum,delt,dt , asum, meanv, mti me, del 1 ,del2, 

: const ,t, sum, sunb, val , sval , wdf (dp*2, dp), coef ,df, AREA 

complex s(dp*2),dun,c(dp*5),duTi3,dLin2 
character*25 inname 

common/ in/ mvopt , redopt ,dimt ,dimf ,mm,nn, rf , rdp, rdp2 # 

: nf ,mt ,nf2,mt2 

common/rl/ pi , tin,ain,rawa,rwdf , rswdf , 

: rawt, dt sun, del t,dt, asum, meanv, mti me, dell ,del2, 

: const, t, sun, sumb,va l , sval, wdf , coef # df # area 

common/comp/ s,dum,c,dum3,dum2 
dp2 = dp* 2 



as Lin = 0.0 

do 100 i = 1 , dp 

asum = a sum + ain(i) 

100 continue 

meanv = asum / dp 

if (mvopt .eq.1 ) then 

do 200 i = 1 , dp 
ain(i) = ain(i) - meanv 
200 continue 

endif 

END 
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SUBROUTINE REDUCE3 



* 



* * 

* This subroutine reduces the wigner ville * 

* to a workable size. * 

* ★ 
******************************************************* 



integer dp,dp2,mvopt ,redopt,dimt,dimf ,rrm,m,rf , rdp,rdp2, 

: nf ,mt ,nf2,mt2 

PARAMETER (DP = 512) 

real pi ,t in(dp),ain(dp),rawa(dp),rwdf(dp*2,dp),rswdf (dp*2,dp), 
: rawt(dp),dtsun,del t,dt,asi*n,meanv,mtime,del 1 ,del2, 

: const, t,sum,sunb,val ,sval ,wdf (dp*2,dp),coef ,df ,area 

complex s(dp*2) ,dum,c(dp*5) ,ckjri5,dum2 
character*25 inname 

common/ in/ mvopt,redopt,dimt,dimf ,mm,nn,rf , rdp, rdp2, 

: nf ,mt,nf2,mt2 

common/rl/ pi ,tin,ain,rawa,rwdf ,rswdf , 

: rawt,dtsum,del t,dt ,asun,meanv,mt i me, dell ,del2, 

: const,t,sun,sunb,val,sval ,wdf ,coef ,df ,area 

common/comp/ s,duri,c,dun3,durn2 
dp2 = dp* 2 



IF (REDOPT.EQ. 1 ) THEN 
RF = 64 
RT = 32 

ELSEIF (REDOPT.EQ. 2) THEN 
RF = 128 
RT = 64 

ELSEIF (REDOPT.EQ. 3) THEN 
RF = 128 
RT = 128 
ELSE 

RF = 256 
RT = 128 

ENDIF 

nn = dp2 / RF 
mm = dp / RT 

l = 0 

do 100 j = 1 , dp2 , nn 
1 = 1 + 1 
k = 0 

do 100 i = 1 , dp , mm 
k = k ♦ 1 

rwdf(l,k) = wdf(j,i) 
100 continue 

END 
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SUBROUTINE SHOOT H3 



************************************************************* 
* ' * 

* This sub rout ir>e reciuces and smooths the WDF along both * 

* the time and frequency axes for purposes of plotting * 

* clarity. * 

* * 

************************************************************* 

* 

real hg(-100: 100, -100:100) 



integer dp,dp2,mvopt , redopt ,dimt ,dimf ,mm,m,rf , rdp.rdpZ, 

: nf ,mt ,nf2,mt2 

PARAMETER (DP = 512) 

real pi ,tin(dp),ain(dp),rawa(dp),rwdf (dp*2,dp),rswdf (dp*2,dp), 
: rawt(dp) ,dtsun,del t ,dt ,asun,meanv,mt ime,del 1 ,del2, 

: const , t,sun,sifTib,val ,sval ,wdf (dp*2,dp) ,coef ,df ,area 

complex s(dp*2),dum,c(dp*5),d<jTi3,dcin2 
character*25 inname 

common/ in/ mvopt, redopt, d i mt, d imf , mm, nn,rf,rdp,rdp2, 

: nf ,mt,nf2,mt2 

common/rl/ pi , tin,ain, rawa,rwdf ,rswdf , 

: rawt,dtsun,delt,dt,asun,meanv,mtime,del1 ,del2, 

: const, t,sun,simb,val ,sval,wdf ,coef ,df ,area 

common/comp/ s,dun,c,dum3,dum2 
dp2 = dp* 2 



df =1 ./(4.*dp*dt) 

rdp=dp/mm 

rdp2=dp2/nn 

nf2=nf*2 

mt2=mt*2 

val=1 ./( (2.*pi )**2*nf*mt) 

do 20 j=-mt2,mt2 
do 10 i=-nf2,nf2 

coef=-(( j*j )/(2.*mt*mt)+( i*i )/(2.*nf*nf )) 
hg( i , j )=val*exp(coef ) 

1 0 cont i nue 
20 cont i nue 

do 100 j=1,rdp 
do 100 i=1,rdp2 
100 rswdf(i, j)=0.0 

* 

* 

i i =0 

DO 500 N=1 ,DP2,nn 

jj = 0 

ii = ii+1 

DO 450 M = 1 ,DP,mm 

jj = jj+1 

DO 400 K=m-MT,rrHMT 

DO 350 L = n-nf,r>+nf 

IF (L.LT.1) THEN 

RSWDFC i i , j j ) = RSWDF(iiJj) ♦ 0 
ELSEIF (K.LT.1) THEN 
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RSUDF(iiJj) = RSWDF(iiJj) + 0 
ELSEIF (L.GT.DP2) THEN 

RSWD F(l l,jj) = RSUOF( l i , j j ) ♦ 0 
ELSEIF (K.GT.DP) THEN 

RSWDFC l i , j j > = RSUDF(n,jj) ♦ 0 
ELSE 

RSUDF(iiJj) = RSUOF( 1 1 , j j )+W)F(L,K)*HG(L-N,K-M) 
ENDIF 



350 

400 



CONTINUE 

CONTINUE 

RSUOF(ii,jj)=RSUDF(ii,jj) 



450 

500 



CONTINUE 

CONTINUE 



END 
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SUBROUTINE WIGH3 



* 

********************************************************* 
* * 

* This subroutine calculates the VCF of the signal * 

* * 
********************************************************* 
* 



integer dp,dp2,mvopt ,redopt,dimt,dimf ,mm,nn, rf ,rdp,rdp2, 

: nf ,mt,nf2,mt2 

PARAMETER (DP = 512) 

real pi , t i n(dp) , ain(dp) , rawa(dp) , rwdf (dp*2 ,dp) , rswdf (dp*2 ,dp) , 
: rawt(dp),dtsLm,delt,dt,asun,meanv,mtime,del1,del2, 

: const, t,sun,sunb,val ,sval ,wdf (dp*2,dp),coef ,df ,area 

complex s(dp*2),dun,c(dp*5),dum3,dun2 
character*25 inname 

common/ in/ mvopt, redopt,dimt,dimf ,mm,nn, rf ,rdp,rdp2, 

: nf ,mt,nf2,mt2 

common/rl/ pi , tin, ain, rawa, rwdf , rswdf , 

: rawt,dtsijn,delt,dt,asum,meanv,mtime,del1 ,del2, 

: const, t ,sun,sumb,val ,sval ,wdf ,coef ,df ,area 

common/comp/ s,dun,c,dtm3,duTi2 
dp2 = dp*2 



do 100 j = 1 , dp 

*******auto correlation of the signal corr3 
coef = 2.0 * dt 
do 300 i = 1 , dp+1 

if(j.ge.i) then 
dum = s( j-i+1) 
else 

dum = cmplx(0.,0.) 
endif 

c(i) = coef * (s( j+i -1 )*conjg(dum)) 

if (i .ne.j.and. i .ne.dp+1 ) then 
c(dp2- i+2) = conjg(c( i )) 
endi f 

300 continue 

call f f t3 

do 200 i=1 ,dp2 

wdf (i , j )=real (c( i ) ) 

200 continue 

1 00 cont i nue 

END 
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C. PLOTTING PROGRAMS 



Listed below are the two computer programs that were used to produce the Figures in 



this work. 

1. 2-D Plotting Routine 



PROGRAM SIGNALPLT 

* 

********************************************************* 



* THIS PROGRAM USES THE GRAPHIC PACKAGE DISSPLA TO 

* PLOT THE A 2-D FUNCTION. 

* 

********************************************************* 



* ***.. -DECLARING VARIABLES---*** 

* 

REAL XARAY( 1000) # YARAY( 1000) 

REAL X, Y # F,XSTART,XEND, YSTART,YEND,XSTEP 
INTEGER I PROPT , NPNTS, IX, I Y 

* 

* *** — INITIALIZING VARIABLES AND SETTING GRAPH OPTIONS---*** 

* 

* — set option for x starting & ending value — 
URITEC*,*), 'DECIDE THE STARTING AND ENDING VALUES OF X' 
URITEC*,*),' INPUT THE STARTING VALUE OF X' 

READC*,*), XSTART 

URITE(*,*), ‘INPUT THE ENDING VALUE OF X' 

READC*,*), XEND 

* — set desired x step size — 

URITEC*,*), 'DECIDE THE X STEP SIZE' 

URITEC*,*), 1 INPUT THE X STEP SIZE' 

READC*,*), XSTEP 

* — set option for y starting & ending value — 
URITEC*,*), 'DECIDE THE STARTING AND ENDING VALUES OF Y' 
URITEC*,*), 'INPUT THE STARTING VALUE OF Y' 

READC*,*), YSTART 

URITEC*,*), 'INPUT THE ENDING VALUE OF Y' 

READC*,*), YEND 

* — set option for number of points to plot — 
URITEC*,*), 'DECIDE THE NUMBER OF POINTS TO PLOT' 
URITEC*,*),' INPUT THE NUMBER OF POINTS' 

READC*,*), NPNTS 

* — set option for line style ( I MARK) — 

URITEC*,*), 'DECIDE THE IMARK STYLE DESIRED' 

URITEC*,*), 'INPUT THE IMARK STYLE' 

READC*,*), IMARK 

* — set option for grid style — 

URITEC*,*), 'DECIDE THE GRID STYLE DESIRED' 

URITEC*,*), 'INPUT THE X GRID STYLE' 

READC*,*), IX 

URITEC*,*), 'INPUT THE Y GRID STYLE' 

READC*,*), IY 

* 

* ---set option for output device--- 

URITEC*,*), 'YOU UILL NOU DECIDE UHERE TO SEND THE OUTPUT' 
URITEC*,*), 'ENTER 0 FOR OUTPUT TO THE SCREEN' 

URITEC*,*), 'ENTER 1 FOR OUTPUT TO THE LASER' 

READC*,*), IPROPT 
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C CALL PDEVC 1 LN03 1 , I EER) 

IF (IPROPT.EQ.O) THEN 
CALL PGPX 
ELSE 

CALL LN03I 
END IF 



* "FUNCTION EXECUTION TO DETERMINE PLOTTING POINTS** 

* 

OPEN (20, F I LE= ' S I N400256A2.DAT ' , STATUS= ' OLD ' ) 

CALL F(XARAY,YARAY,NPNT$,X$TART,XEND) 



* ***-- -LEVEL ONE COMMANDS---*** 

CALL HWROT( 'COMIC') 

CALL PAGE(1 1 .0,8.5) 

CALL NOBRDR 

CALL AREA2D(9. 0,6.0) 

* 

* ***--- LEVEL TWO COMMANDS---*** 

* 

CALL SUISSB 

CALL HEADINC400 HZ SINE WAVES* , 100, 1 .5, 2) 

CALL HEADINCSIGNAL LENGTH = 0.128 SEC / AMP = 2. OS* , 100, 1 .5, 2) 
CALL SWISSM 

CALL XNAME( 1 TIME (sec)S',100) 

CALL YNAMEC ‘AMPLITUDES' ,100) 

CALL GRAF(XSTART,XSTEP,XEND, YSTART, 'SCALE * ,YEND) 

* 

* ***--- LEVEL THREE COMMANDS---*** 

* 

CALL CURVE ( XARAY, YARAY,NPNTS, IMARK) 

CALL FRAME 
CALL GRID( IX, I Y) 

* 

* *** CLOSING COMMANDS---*** 

* 

CALL ENDPL(O) 

CALL DONEPL 
END 



* ---FUNCTION SUBROUTINE TO DETERMINE PLOTTING POINTS--- 

★ 

SUBROUTINE F ( XARAY, YARAY, NPNTS,XSTART,XEND) 

* 

REAL XARAY(NPNTS) , YARAY (NPNTS) ,X$TART,XEND, STEP, XTEMP 
INTEGER NPNTS 

* ---READ IN THE XARAY AND YARAY VALUES--- 

DO 100 I = 1, NPNTS 

READ(20,*) XARAY(I), YARAY (I ) 

100 CONTINUE 
END 
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2. 3-D Plotting Routine 



PROGRAM UDFPLTVER3 

* 

********************************************************* 



* THIS PROGRAM USES THE GRAPHIC PACKAGE DISSPLA TO 

* PLOT A 3-D FUNCTION. 

* 

********************************************************* 



* *** DECLARING VARIABLES---*** 

* 

REAL RWDF<32768) 

INTEGER IPROPT, IXPTS, IYPTS 

* 

* *** INITIALIZING VARIABLES AND SETTING GRAPH OPTIONS---*** 

OPEN(15,FILE='RSWDF.OUT',STATUS='OLD') 

* 

DO 100 I = 1,32768 
READ(15,*) RWDF(I) 

100 CONTINUE 



call pdev('ln03', ieer) 

* 

* **FUNCT ION EXECUTION TO DETERMINE PLOT** 

* 

* ***-.- LEVEL ONE COMMANDS---*** 

* 

CALL PAGE(8.5, 11.0) 

CALL AREA2D(8. 0,7.0) 

* 

* *** — LE VEL TWO COMMANDS---*** 

* 

CALL VOLM3D(6.,6.,4.) 

CALL HEADINC200 & 600 Hz SINE WAVES' , 100, -1 .5,2) 

CALL HEAD IN( 'FREQUENCY RESOLUTION ; 5 X 35S • , 100, - 1 .5 ,2) 
CALL X3NAME( ' FREQUENCYS ' , 100) 

CALL Y3NAMEC ',1) 

CALL Z3NAMEC 'AMPLITUDES', 100) 

CALL VUANGL( -90. ,0. ,30. ) 

CALL GRAF3D(0.,1250./5., 1250.0, 

: 0.,. 2048/2. ,0.2048, 

: 0., 'SCALE', .01) 

* 

* ***... level three commands---*** 

* 

CALL SURMAT(RWDF, 1,256, 1,128,0) 

CALL FRAME 

* 

* *** CLOSING COMMANDS---*** 

* 

CALL ENDPL(O) 

CALL DONEPL 

* 

END 
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